% input contains
%     numg   = number of energy groups
%     numm   = number of materials
%     xcm    = coarse divisions
%     xfm    = fine mesh interval for coarse divisions
%     mt     = material assignment for each coarse division
%     data   = cross-sections in the form
%              (mat1/g1) sigTOT   sigA  sigSg1->g1  sigSg2->g1 ...
%                 ...g2) sigTOT   sigA  sigSg1->g2  ....
%              (mat2/g1) ...
%     src    = volumetric (isotropic) source by coarse mesh and energy
%              (for adj, src is the cross-section for the source region)
%     ord    = number of ordinates (2,4,8, or 12)
%     maxit  = maximum iterations
%     maxerr = maximum relative pointwise error in phi

format short

xcm    = [ 0  1.26];
xfm    = [  60 ];
ycm    = [ 0  1.26];
yfm    = [  60 ];
src(1,:,:)    = 1;            
mt     = 1;
           % St   Sa  
data   = [   1.0  1.  0.0];


input   =   struct(   ...
    'numg',            1, ...     % number of groups
    'numm',            1, ...     % number of materials
    'xcm',           xcm, ...     % slab bounds
    'xfm',           xfm, ...     % number of fine meshes
    'ycm',           ycm, ...     % slab bounds
    'yfm',           yfm, ...     % number of fine meshes    
    'mt',             mt, ...     % slab material ids
    'data',         data, ...     % mat comp's
    'src',           src, ...     % volume source
    'ord',            50, ...     % number of ordinates
    'quad',            2, ...     % quad type; 1=level sym. 2 = UEN
    'maxit',          1, ...     % max iterations
    'maxerr',       1e-8, ...     % max pointwise phi error
    'adj',             0  ...     % adjoint flag
    );

%[phi,psi,psiV,psiH] = sn_two_d(input);
%[KK,psi2,psiH,psiV]=sn_two_d_matrix(input);

%[phi] = sn_two_d(input);
[phi] = sn_two_d_mem(input);

% plot phi on the coarse grid
numx = length(xcm)-1;
numy = length(ycm)-1;
phi_coarse = zeros(numx,numy);
% number of fine meshes in x and y coordinates
N = sum(xfm);
M = sum(yfm);
% number of coarse meshes in x and y coordinates
CN = length(xfm);
CM = length(yfm);
% compute dx and dy vectors
dx = zeros(N,1);
dy = zeros(M,1);
% coarse mesh index for fine mesh
cix=zeros(N,1);
ciy=zeros(M,1);
j = 0;
for i = 1:CN
    dx( (j+1):(j+xfm(i)) ) = (xcm(i+1)-xcm(i))/xfm(i);    
    cix( (j+1):(j+xfm(i)) )=i;
    j = sum(xfm(1:i));
end
j = 0;
for i = 1:CM
    dy( (j+1):(j+yfm(i)) ) = (ycm(i+1)-ycm(i))/yfm(i);    
    ciy( (j+1):(j+yfm(i)) )=i;
    j = sum(yfm(1:i));
end
for j = 1:N
    for i = 1:M
        phi_coarse( cix(i), ciy(j) ) = phi_coarse( cix(i), ciy(j) ) + ...
                                       phi(i,j)*dx(i)*dy(i);
    end
end
for j = 1:numy
    for i = 1:numx
        V(i,j) = (xcm(i+1)-xcm(i))*(ycm(j+1)-ycm(j));
    end
end
phi_coarse = phi_coarse ./ V
    

%surf( psiH(:,:,1,1) );
%spy( abs(K-KK)>0.000001
%surf(phi)